On the Generation of Positivstellensatz Witnesses 
in Degenerate Cases* 

David Monniaux^ Pierre Corbineau^ 
January 22, 2013 

Abstract 

One can reduce the problem of proving that a polynomial is nonneg- 
ative, or more generally of proving that a system of polynomial inequal- 
ities has no solutions, to finding polynomials that are sums of squares 
of polynomials and satisfy some linear equality {Positivstellensatz) . This 
produces a witness for the desired property, from which it is reasonably 
easy to obtain a formal proof of the property suitable for a proof assistant 
such as Coq. 

The problem of finding a witness reduces to a feasibility problem in 
semidefinite programming, for which there exist numerical solvers. Un- 
fortunately, this problem is in general not strictly feasible, meaning the 
solution can be a convex set with empty interior, in which case the nu- 
merical optimization method fails. Previously published methods thus 
assumed strict feasibility; we propose a workaround for this difficulty. 

We implemented our method and illustrate its use with examples, 
including extractions of proofs to Coq. 

1 Introduction 

Consider the following problem: given a conjunction of polynomial equalities, 
and (wide and strict) polynomial inequalities, with integer or rational coeffi- 
cients, decide whether this conjunction is satisfiable over M; that is, whether 
one can assign real values to the variables so that the conjunction holds. A 
particular case is showing that a given polynomial is nonnegative. 

The decision problem for real polynomial inequalities can be reduced to 
quantifier elimination: given a formula F, whose atomic formulas are polyno- 
mial (in)equalities, containing quantifiers, provide another, equivalent, formula 
F' , whose atomic formulas are still polynomial (in)equalities, containing no 
quantifier. An algorithm for quantifier elimination over the theory of real closed 
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fields (roughly speaking, (M, 0, 1, +, x, <) was first proposed by Tarski 29l 26|. 
but this algorithm had non-elementary complexity and thus was impractical. 
Later, the cylindrical algebraic decomposition (CAD) algorithm was proposed 
0, with a doubly exponential complexity, but despite improvements f§] CAD 
is still slow in practice and there are few implementations available. 

Quantifier elimination is not the only decision method. Basu et al. 0, The- 
orem 3] proposed a satisfiability testing algorithm with complexity s'^^-'^d'^^'^^ 
where s is the number of distinct polynomials appearing in the formula, d is their 
maximal degree, and k is the number of variables. We know of no implemen- 



tation of that algorithm. Tiwari [31[ proposed an algorithm based on rewriting 
systems that is supposed to answer in reasonable time when a conjunction of 
polynomial inequalities has no solution. 

Many of the algebraic algorithms are complex, which leads to complex im- 
plementations. This poses a methodology problem: can one trust their results? 
The use of computer programs for proving lemmas used in mathematical theo- 
rems was criticized in the case of Thomas Hales' proof of the Kepler conjecture. 
Similarly, the use of complex decision procedures (as in the proof assistant 
PV^ or program analyzers (as, for instance, AstrecQ) in order to prove the 
correctness of critical computer programs is criticized on grounds that these 
verification systems could themselves contain bugs. 

One could formally prove correct the implementation of the decision proce- 
dure using a proof assistant such as Coq 3^, 1^ ; but this is likely to be long and 



difficult. An alternative is to arrange for the procedure to provide a witness of 
its result. The answer of the procedure is correct if the witness is correct, and 
correctness of the witness can be checked by a much simpler procedure, which 
can be proved correct much more easily. 

Unsatisfiability witnesses for systems of complex equalities or linear rational 
inequalities are already used within DPLL(r) satisfiability modulo theory deci- 
sion procedures [lil . ch. Il][l3|. It is therefore tempting to seek unsatisfiability 
witnesses for systems of polynomial inequalities. 

In recent years, it was suggested [13 to use numerical semidefinite program- 
ming to look for proof witnesses whose existence is guaranteed by a Positivstel- 



lensatz [15|, |28j, |25| . The original problem of proving that a system of polynomial 
inequalities has no solution is reduced to: given polynomials Pi and R, derived 
from those in the original inequalities, find polynomials Qi that are sums of 
squares such that Yli PiQi = R- Assuming some bounds on the degrees of Qi, 
this problem is in turn reduced to a semidefinite programming pure feasibility 



problem [6|, |32| , a form of convex optimization. The polynomials Qi then form 
a witness, from which a machine- checkable formal proof, suitable for tools such 
as Coq "SO^ or Isabelle flTl, may be constructed. 

Unfortunately, this method suffers from a caveat: it applies only under a 
strict feasibility condition [2l|: a certain convex geometrical object should not 
be degenerate, that is, it should have nonempty interior. Unfortunately it is 
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very easy to obtain problems where this condition is not true. Equivalently, the 
method of rationahzation of certificates 13[ has a limiting requirement that the 
rationalized moment matrix remains positive semidefinite. 

In this article, we explain how to work around the degeneracy problem: we 
propose a method to look for rational solutions to a general SDP feasibility 
problem. We have implemented our method and applied it to some examples 
from the literature on positive polynomials, and to examples that previously 
published techniques failed to process. 



2 Witnesses 

For many interesting theories, it is trivial to check that a given valuation of the 
variables satisfies a quantifier-free formula. A satisfiability decision procedure 
will in this case tend to seek a satisfiability witness and provide it to the user 
when giving a positive answer. 

In contrast, if the answer is that the problem is not satisfiable, the user has 
to trust the output of the satisfiability testing algorithm, the informal meaning 
of which is "I looked carefully everywhere and did not find a solution." In 
some cases, it is possible to provide unsatisfiability witnesses: solutions to some 
form of dual or auxiliary problem that show that the original problem had no 
solution. 



2.1 Nonnegativity Witnesses 

To prove that a polynomial P is nonnegative, one simple method is to express 
it as a sum of squares of polynomials. One good point is that the degree of the 
polynomials involved in this sum of squares can be bounded, and even that the 
choice of possible monomials is constrained by the Newton polytope of P, as 
seen in §31 

Yet, there exist nonnegative polynomials that cannot be expressed as sums 
of squares, for instance this example due to Motzkin [23| : 

A/f— 6,42,24 q222 

Ivl — ' ^2*^3 ' "^2^3 — yi- ) 

However, Artin's answer to Hilbert's seventeenth problem is that any non- 
negative polynomial can be expressed as a sum of squares of rational functions^ 

It follows that such a polynomial can always be expressed as the quotient 
Q2/Q1 of two sums of squares of polynomials, which forms the nonnegativity 
witness, and can be obtained by solving P.Qi — Q2 — ioi Qi (this result 
is also a corollary of Th. [1]). 

^ There exists a theoretical exact algorithm for computing such a decomposition for homo- 
geneous polynomials of at most 3 variables [l4j : we know of no implementation of it and no 
result about its practical usability. 
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2.2 Unsatisfiability Witnesses for Polynomial Inequalities 

For the sake of simplicity, we shall restrict ourselves to wide inequalities (the 
extension to mixed wide/strict inequalities is possible). Let us first remark that 
the problem of testing whether a set of wide inequalities with coefficients in a 
subfield K of the real numbers is satisfiable over the real numbers is equivalent to 
the problem of testing whether a set of equalities with coefficients K is satisfiable 
over the real numbers: for each inequality P{xi, . . . ^Xm) > 0, replace it by 
P{xi, . . . , Xm) — — ^, where /j, is a new variable. Strict inequalities can also 
be simulated as follows: Pi{xi, . . . , Xm) 7^ is replaced by Pi{xi, . . . , Xm)-fJ- = 1 
where is a new variable. One therefore does not gain theoretical simplicity by 
restricting oneself to equalities. 



Stengle [28| proved two theorems regarding the solution sets of systems of 
polynomial equalities and inequalities over the reals (or, more generally, over 
real closed fields): a Nullstellensatz and a Positivstellensatz; a similar result 
was proved by Krivine [l5j. Without going into overly complex notations, let us 
state consequences of these theorems. 

Let K be an ordered field (such as Q) and K' be a real closed field containing 
K (such as the real field R), and let X be a list of variables Xi,. . . A*"^ 
denotes the squares of elements of A. The multiplicative monoid generated by 
A is the set of products of zero of more elements from A. The ideal generated 
by A is the set of sums of products of the form PQ where Q G K[K.] and P G A. 
The positive cone generated by A is the set of sums of products of the form 
p.P.Q^ where p G K , p > 0, P is in the multiplicative monoid generated by A, 
and Q € KlX]. Remark that we can restrict P to be in the set of products of 
elements of A where no element is taken twice, with no loss of generality. 

The result 18, 17, ^ of interest to us is: 



Theorem 1. Let Fy, F>, F^ be sets of polynomials in ^^[X], to which we 
impose respective sign conditions > 0, > 0, = 0, 0. The resulting system is 
unsatisfiable over K'"' if and only if there exist an equality in K[K] of the type 
S + P + Z = 0, with S in the multiplicative monoid generated by P> U F^ , P 
belongs to the positive cone generated by F> U Fy , and Z belongs to the ideal 
generated by F^ . 

{S, P, Z) then constitute a witness of the unsatisfiability of the system. 
For a simple example, consider the following system, which obviously has no 
solution: 

-2 + y2 > 
1 - > 



(2) 



A Positivstellensatz witness is y^{—2 + y^) + 1(1 — y^) + + 1 = 0. Another 
is (I + t) {~2 + y^) + ^{l-y^) + l = Q. 

Consider the conjunction C: Pi > OA- ■ -APn > where Pi £ Q[Xi, . . . , X„i]- 
Consider the set H of products of the form J|,- P^"^ for w G {0, 1}" — that is, 



^Another result, due to Schmiidgen |25|1 . gives simpler witnesses for Pi > A • ■ ■ A Pn > 
=> C in the case where Pi > A ■ ■ ■ A P„ > defines a compact set. 
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the set of all products of the Pi where each Pi appears at most once. Obviously, 
if one can exhibit nonnegative functions Qj such that X^r en Qj^^j + 1 = 0, then 
C does not have solutions. Theorem [T] guarantees that if C has no solutions, 
then such functions Qj exist as sum of squares of polynomials (we simply apply 
the theorem with = = ^ and thus S = {!}). We have again reduced our 
problem to the following problem: given polynomials Tj and R, find sums-of- 
squares polynomials Qj such that QjTj = R. Because of the high cost of 
enumerating all products of the form Y\i Pj"^ , we have first looked for witnesses 
of the form J2TjeS Qj^J + 1 = 0- 

3 Solving the Sums-of- Squares Problem 

In i i2.1l and W2.21 we have reduced our problems to: given polynomials {Pj)i<j<n 
and R in Q[Xi, . . . , X.^], find polynomials that are sums of squares Qj such that 

J2pjQj = r (3) 



We wish to output the Qj as Qj = X]r=i where aji G Q"*" and Lji are 

polynomials over Q. We now show how to solve this equation. 

3.1 Reduction to Semidefinite Programming 

Lemma 1. Let P G Y, . . .] be a sum of squares of polynomials P^ . Let 

M = {mi, . . . , TO|jvfj} be a set such that each Pi can be written as a linear combi- 
nation of elements of M (M can be for instance the set of monomials in the Pi). 
Then there exists a \M\ x \M\ symmetric positive semidefinite matrix Q with 
coefficients in K such that P{X, Y, . . . ) = [mi, . . . , m|jvf • • ■ : ^\m\]'^ , not- 

ing the transpose of v. 

Assume that we know the ALj, but we do not know the matrices Qj. The 
equality Pj{MjQj{Mj)'^) — R directly translates into a system [S) of affine 

linear equalities over the coefficients of the Qj: J2ji^'^jQji^j)'^)-Pj — ^ is the 
zero polynomial, so its coefficients, which are affine linear combinations of the 
coefficients of the Qj matrices, should be zero; each of these combinations thus 
yields an affine linear equation. The additional requirement is that the Qj are 
positive semidefinite. 

One can equivalently express the problem by grouping these matrices into a 
block diagonal matrix Q and express the system (S) of affine linear equations 
over the coefficients of Q. By exact rational linear arithmetic, we can obtain a 
system of generators for the solution set of (S): Q £ —Fq + vect(Fi, . . . , Fm). 
The problem is then to find a positive semidefinite matrix within this search 
space; that is, find ai, . . . , a„i such that —Fo + J2i '^i^i — 0- This is the problem 
of semidefinite programming: finding a positive semidefinite matrix within an 
affine linear variety of symmetric matrices, optionally optimizing a linear form 
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For instance, the second unsatisfiability witness we gave for constraint sys- 
teni[2]is defined, using monomials {1,?/}, 1 and {f,y}, by: 




T 



2. 



\ 




0/ 



It looks like finding a solution to Equ. [3] just amounts to a SDP problem. 
There are, however, several problems to this approach: 

1. For the general Positivstellensatz witness problem, the set of polynomials 
to consider is exponential in the number of inequalities. 

2. Except for the simple problem of proving that a given polynomial is a sum 
of squares, we do not know the degree of the Qj in advance, so we cannot 
choose finite sets of monomials Mj . The dimension of the vector space 
for Qj grows quadratically in |Afj|. 

3. Some SDP algorithms can fail to converge if the problem is not strictly 
feasible — that is, the solution set has empty interior, or, equivalently, is 
not full dimensional (that is, it is included within a strict subspace of the 
search space). 

4. SDP algorithms are implemented in floating-point. If the solution space is 
not full dimensional, they tend to provide solutions Q that are "almost" 
positive semidefinite (all eigenvalues greater than — e for some small posi- 
tive e), but not positive semidefinite. 

Regarding problem 1, bounds on degrees only matter for the completeness 
of the refutation method: we are guaranteed to find the certificate if we look in 
a large enough space. They are not needed for soundness: if we find a correct 
certificate by looking in a portion of the huge search space, then that certificate 
is correct regardless. This means that we can limit the choice of monomials 
in Mj and hope for the best. 

Regarding the second and third problems : what is needed is a way to reduce 
the dimension of the search space, ideally up to the point that the solution set 
is full dimensional. As recalled by [2lj, in a sum-of-square decomposition of a 
polynomial P, only monomials x"^ . . . x"" such that 2(ai, . . . , a„) lies within the 
Newton polytop^of P can appear [22|, Th. 1]. This helps reduce the dimension if 
P is known in advance (as in a sum-of-squares decomposition to prove positivity) 
but does not help for more general equations. 

Kaltofen et al. suggest solving the SDP problem numerically and looking 
for rows with very small values, which indicate useless monomials that can be 

^ There exist non-elementary bounds on the degree of the monomials needed [T3l . In the 
case of Schmiidgen's result on compact sets, there are better bounds [25j. 

^The Newton polytope of a polynomial P, or in Reznick's terminology, its cage, is the 
convex hull of the vertices (oi, . . . , a„) such that ^ . . monomial of P. 
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safely removed from the basis; in other words, they detect "approximate kernel 
vectors" from the canonical basis. Our method is somehow a generalization 
of theirs: we detect kernel vectors whether or not they are from the canonical 
basis. 

In the next section, we shall investigate the fourth problem: how to deal 
with solution sets with empty interior. 

3.2 How to Deal with Degenerate Cases 

In the preceding section, we have shown how to reduce the problem of finding 
unsatisfiability witnesses to a SDP feasibility problem, but pointed out one 
crucial difficulty: the possible degeneracy of the solution set. In this section, we 
explain more about this difhculty and how to work around it. 

Let K. be the cone of positive semidefinite matrices. We denote by M ^ 
a positive semidefinite matrix Af , by Af >- a positive definite matrix M . The 
vector y is decomposed into its coordinates yi. x denotes a floating-point value 
close to an ideal real value x. 

We consider a SDP feasibility problem: given a family of symmetric matrices 
FQ,Fi, . . . ,Fm, find {yi)i<i<m such that 



The Fi have rational coefficients, and we suppose that there is at least one 
rational solution for y such that F{y) >z_ 0. The problem is how to find such a 
solution. 

If nonempty, the solution set S C R™ for the y, also known as the spectra- 
hedron, is semialgebraic, convex and closed; its boundary consists in y defining 
singular positive semidefinite matrices, its interior are positive definite matri- 
ces. We say that the problem is strictly feasible if the solution set has nonempty 
interior. Equivalently, this means that the convex S has dimension m. 

Interior point methods used for semidefinite feasibility, when the solution set 
has nonempty interior, tend to find a solution y in the interior away from the 
boundary. Mathematically speaking, if y is a numerical solution in the interior 
of the solution set, then there is e > such that for any y such that ||y — y|| < e, 
y is also a solution. Choose a very close rational approximation y of y, then 
unless we are unlucky (the problem is almost degenerate and all any suitable e 
is extremely small), then y is also in the interior of S. Thus, F{y) is a solution 
of problem |4l 

This is why earlier works on sunis-of-square methods [2l| have proposed 
finding rational solutions only when the SDP problem is strictly feasible. In 
this article, we explain how to do away with the strict feasibility clause. 

Some problems are not strictly feasible. Geometrically, this means that the 
hnear affine space {— Fo -I- J2^i Vi^i I (j/Ij ■ ■ • '2/™) ^ K™} is tangent to the 
semidefinite positive cone /C. Alternatively, this means that the solution set is 
included in a strict linear affine subspace of R™ . Intuitively, this means that we 




(4) 
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are searching for the solution in "too large a space" ; for instance, if m = 2 and 
y lies in a plane, this happens if the solution set is a point or a segment of a line. 
In this case, some SDP algorithms may fail to converge if the problem is not 
strictly feasible, and those that converge, in general, will find a point slightly 
outside the solution set. The main contribution of this article is a workaround 
for this problem. 

3.3 Simplified algorithm 

We shall thus now suppose the problem has empty interior. 
The following result is crucial but easily proved: 

Lemma 2. Let E be a linear affine subspace of the n x n symmetric matrices 
such that E f] K. ^ %. F in the relative interior I of E (1 IC. Then it follows: 

1. For all F' e Er\IC,kerF CkerF'. 

2. The least affine space containing EC\K, is H ~ {AI G E \ ker Af D keri^}. 

Suppose we have found a numerical solution y, but it is nearly singular — 
meaning that it has some negative eigenvalues extremely close to zero. This 
means there is v ^ such that |v.F(y)| < e||v||. Suppose that y is very 
close to a rational solution y and, that v.F{y) = 0, and also that y is in the 
relative interior of 5' — that is, the interior of that set relative to the least 
linear affine space containing S. Then, by lemma El all solutions F{y') also 
satisfy v.F(y') = 0. Remark that the same lemma implies that either there is 
no rational solution in the relative interior, or that rational solutions are dense 
in S. 

How can finding such a v help us? Obviously, if v G 01=0 ^'^^ ^» ' discovery 
does not provide any more information than already present in the linear affine 
system — Fq + Vect {Fi, . . . , F„j). We thus need to look for a vector outside that 
intersection of kernels; then, knowing such a vector will enable us to reduce the 
dimension of the search space from m to m' < m. 

Thus, we look for such a vector in the orthogonal complement of fXiLo ker i^j, 
which is the vector space generated by the rows of the symmetric matrices 
Fq, . . . , Fm ■ We therefore compute a full rank matrix B whose rows span the 
exact same space; this can be achieved by echelonizing a matrix obtained by 
stacking Fq, . . . , F„i. Then, v — wB for some vector w. We thus look for w 
such that G(y).w = 0, with G(y) = BF{y)B'^. 

The question is how to find such a w with rational or, equivalently, inte- 
ger coefficients. Another issue is that this vector should be "reasonable" — it 
should not involve extremely large coefficients, which would basically amplify 
the floating-point inaccuracies. 

We can reformulate the problem as: find w e Z'" \ {0} such that both w 
and G(y).w are "small", two constraints which can be combined into a single 
objective to be minimized Q;^||G(y).w||2 + ||w|j2, where a > is a coefficient 
for tuning how much we penalize large values of ||G(y).w||2 in comparison to 
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large values of ||w||2. If a is large enough, the difference between aG{y) and 
its integer rounding M is small. We currently choose a — ao/||G(y)||, with 
||M|| the Frobenius norm of M (the Euclidean norm for n x n matrices being 
considered as vectors in R" ), and ao = 10^^. 

We therefore try searching for a small (with respect to the Euclidean norm) 
nonzero vector that is an integer linear combination of the li — (0, . . . , 1, . . . , 0, m^) 
where rrii is the i-th row of M and the 1 is at the i-th position. Note that, be- 
cause of the diagonal of ones, the li form a free family. 

This problem is known as finding a short vector in an integer lattice, and 
can be solved by the Lenstra-Lenstra-Lovasz (LLL) algorithm. This algorithm 
outputs a free family of vectors Si such that si is very short. Other vectors in 
the family may also be very short. 

Once we have such a small vector w, using exact rational linear algebra, we 
can compute Fq, . . . , F^^, such that 



The resulting system has lower search space dimension m' < in, yet the same 
solution set dimension. By iterating the method, we eventually reach a search 
space dimension equal to the dimension of the solution set. 

If we find no solution Fq, then it means that the original problem had no 
solution (the Positivstcllcnsatz problem has no solution, or the monomial bases 
were too small), or that a bad vector v was chosen due to lack of numerical 
precision. This is the only bad possible outcome of our algorithm: it may fail 
to find a solution that actually exists; in our experience, this happens only on 
larger problems (search space of dimension 3000 and more), where the result is 
sensitive to numerical roundoff. In contrast, our algorithm may never provide a 
wrong result, since it checks for correctness in a final phase. 

3.4 More Efficient Algorithm 

In lieu of performing numerical SDP solving on _F = —Fa + '^yiFi ^ 0, we 
can perform it in lower dimension on —{BFqB^) + J^ViiBFiB^) >z 0. Recall 
that the rows of B span the orthogonal complement of CliLQ^ei Fi, which is 
necessarily included in ker F; we are therefore just leaving out dimensions that 
always provide null eigenvalues. 

The reduction of the sums-of-squares problem (Eq. [S]) provides matrices 
with a fixed block structure, one block for each Pj: for a given problem all 
matrices Fq, Fi, . . . , Fm are block diagonal with respect to that structure. We 
therefore perform the test for positive semidefiniteness of the proposed F{y) 
solution block-wise (see Sec. l3.6l for algorithms). For the blocks not found to be 




i^o + ^ 1 (2/1 , . . . , y™) e J n {F I F.v = 0} (5) 
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positive semidefinite, the corresponding blocks of the matrices B and F{y) are 
computed, and LLL is performed. 

As described so far, only a single v kernel vector would be supplied by LLL 
for each block not found to be positive semidefinite. In practice, this tends to 
lead to too many iterations of the main loop: the dimension of the search space 
does not decrease quickly enough. We instead always take the first vector v^^^ of 
the LLL-reduced basis, then accept following vectors v^*^ if ||v(*^||i < /3.||v(^^||i 
and ||G(y).v(') II2 < 7.||G'(y).v(-^) II2. For practical uses, we took /? = 7 = 10. 

When looking for the next iteration y', we use the y from the previous 
iteration as a hint: instead of starting the SDP search from an arbitrary point, 
we start it near the solution found by the previous iteration. We perform least- 
square minimization so that —Fq + X]fc=i Vi-^l the best approximation of 

--Fb + Eili y*^* I (yi,---,2/m)- 



3.5 Extensions and Alternative Implementation 



As seen in ijH our algorithm tends to produce solutions with large numerators 
and denominators in the sum-of-square decomposition. We experimented with 
methods to get F{y') « F{y) such that F{y') has a smaller common denomi- 
nator. This reduces to the following problem: given v S fo -I- vect(fi, . . . , f„) a 
real (floating-point) vector and fp, . . . , f„ rational vectors, find y'j^, . . . , y„ such 
that v' = fo + J2i y'i^i ~ ^ the numerators of v' have a tunable magnitude 
(parameter fi). One can obtain such a result by LLL reduction of the rows of: 



M = 



/Z(/3A*(fo - v)) Z(/3fo) 1 
Z(^^fi) Z(^fi) 1 



V Z(/3Mfn 



Z(/3f„) 



1/ 



(6) 



where /? is a large parameter (say, 10^^) and Z(i;) stands for the integer rounding 
of V. After LLL reduction, one of the short vectors in the basis will be a 
combination Vi^i where lQ,...,ln are the rows of Af, such that i/q ^ 0. 
Because of the large /3// coefficient, yo(fo — v) -I- J27=i Vi^i should be very small, 
thus /o + X]r=i Vi^i ~ ^- -'^^t among those vectors, the algorithm chooses one 
such that Y^^=oyA ^lo* large — and among the suitable v' , the vector of 
numerators is proportional to X]i=o Di^i- 

After computing such a y', we check whether F{y') > 0; we try this for a 
geometrically increasing sequence of /i and stop as soon as we find a solution. 
The matrices Qj then have simpler coefficients than the original ones. Unfor- 
tunately, it does not ensue that the sums of square decompositions of these 
matrices have small coefficients. 

An alternative to finding some kernel vectors of a single matrix would be to 
compute several floating-point matrices, for instance obtained by SDP solving 
with optimization in multiple directions, and find common kernel vectors using 
LLL. 
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3.6 Sub-algorithms and Implementation 

The reduction from the problem expressed in Eq. [3] to SDP with rational solu- 
tions was implemented in SageQ 

Solving the systems of linear equations {S) (Sec. 13.11 over the coefficients 
of the matrices) and [Sj in order to obtain a system —Fq + vect{Fi, . . . ,F,„) 
of generators of the solution space, is done by echelonizing the equation sys- 
tem (in homogeneous form) in exact arithmetic, then reading the solution off 
the echelon form. The dimension of the system is quadratic in the number of 
monomials (on the problems we experimented with, dimensions up to 7900 were 
found); thus efficient algorithms should be used. In particular, sparse Gaussian 
elimination in rational arithmetic, which we initially experimented, is not effi- 



cient enough; we thus instead use a sparse multi-modular algorithm |27l . ch. 7] 
from LinBojO- Multi-modular methods compute the desired result modulo some 
prime numbers, and then reconstruct the exact rational values. 

One can test whether a symmetric rational matrix Q is positive semidefinite 
by attempting to convert it into its Gaussian decomposition, and fail once one 
detects a negative diagonal element, or a nonzero row with a zero diagonal 
element (Appendix. We however experimented with three other methods 
that perform better: 

• Compute the minimal polynomial of Q using a multi-modular algorithm 
[H. The eigenvalues of Q are its roots; one can test for the presence of 
negative roots using Descartes' rule of signs. Our experiments seem to 
show this is the fastest exact method. 

• Compute the characteristic polynomial of Q using a multimodular algo- 
rithm [H and do as above. Somewhat slower but more efficient than Gaus- 
sian decomposition. 

• Given a basis B of the span of Q, compute the Cholesky decomposition 
of B^QB by a numerical method. This decomposition fails if and only if 
B^QB is not positive definite (up to numerical errors), thus succeeds if 
and only if Q is positive semidefinite (up to numerical errors). 

For efficiency, instead of computing the exact basis B of the span of Q, 
we use B from §3.31 whose span includes the span of Q. The only risk is 
that ker B C ker Q while Q is positive semidefinite, in which case B^QB 
will have nontrivial nuUspace and thus will be rejected by the Cholesky 
decomposition. This is not a problem in our algorithm: it just means 
that the procedure for finding kernel vectors by LLL will find vectors in 
kerQ \ ker_B. 

One problem could be that the Cholesky decomposition will incorrectly 
conclude that B^QB is not positive definite, while it is but has very small 



^Sage is a computer algebra system implemented using the Python programming language, 
available under the GNU GPL from http://www.sag6niath.org 

^LinBox is a library for exact linear arithmetic, used by Sage for certain operations. 
|http : //www . llnaig ■ org/ | 
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positive eigenvalues. In this case, our algorithm may then find kernel 
vectors that are not really kernel vectors, leading to an overconstrained 
system and possibly loss of completeness. We have not encountered such 
cases. 

Another problem could be that a Cholesky decomposition is obtained from 
a matrix not positive semidefinite, due to extremely bad numerical behav- 
ior. At worst, this will lead to rejection of the witness when the allegedly 
semidefinite positive matrices get converted to sums of squares, at the end 
of the algorithm. 

Numerical SDP solving is performed using DSD#| [1, i, communicating 
using text files. LLL reduction is performed by fpLLLj^ Least square projection 
is performed using Lapack's DGELS. 

The implementation is available from the first author's Web page ( http : //bit . ly/f BNLhR| 
and |http://bit. ly/gPXNF8) . 



3.7 Preliminary Reductions 

The more coefficients to find there are, the higher the dimension is, the longer 
computation times grow and the more likely numerical problems become. Thus, 
any cheap technique that reduces the search space is welcome. 

If one looks for witnesses for problems involving only homogeneous poly- 
nomials, then one can look for witnesses built out of a homogeneous basis of 
monomials (this technique is implemented in our positivity checker). 

One could also make use of symmetries inside the problem. For instance, if 
one looks for a nonnegativity witness P = N / D of a polynomial P, and P is 
symmetric (that is, there exists a substitution group E for the variables of P such 
that P.a = P for cr G S), then one may reduce the search to symmetric N and 
D. If P = N/D is a witness, then DP = N thus for any cr, {D.a)P = {N.a) and 
thus (X;„ D.(j)P = {Y.a ^■^)^ thus D' = J2a D-a and N' = Y.a ^■'^ constitute 
a symmetric nonnegativity witness. 



4 Examples 

The following system of inequalities has no solution (neither Redlog nor QepCad 
nor Mathcmatica 5 can prove it; Mathematica 7 can): 

Pi = + + -I- z + 1 > 

P2 = - 2y2 + a; + 2 > P3 = + - z > , . 

P4 = -bx^z'^ ~ SOxyz^ - 125y2z3 + 2x'^y'^ + 2Qxy'^ + SOy'* - 2x^ 
-lOx'^y - 2bxy'^ - Ibz^ - Ax^ - 21xy - 47y'^ -3x-y-8>0 

^DSDP is a sdp tool available from http : //www . mcs . a nl ■ gov/DSDP/] 

^"^fpLLL is a LLL library from Damien Stehle et al., available from 
|http : //perso . ens-lyon. f r/damlen. stehle/ | 
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Qi = 8006878^1 + 29138091^2 + 25619868453870/4003439^3 + 14025608^4 + 14385502^^ 
+ 85108577038951965167/12809934226935Ag 

Q2 = 8006878BJ + 25616453B| + 108749058736871/4003439S^ + 161490847987681 
/256I6453B4 + 7272614B| + 37419351S| + 13078817768190/3636307B^ + 71344030945385471151 
/15535579819553B| + 539969700325922707586/161490847987681Bg + 41728880843834 
/12473117BJ(, + 131008857208463018914/62593321265751Bi 1^ , where 

Ai = -1147341/4003439a;j2;3 - 318460/4003439a;^a;3 + xl A2 = X2xl A3 = -4216114037644 
/12809934226935a:ia;3 + 2:23:3 Ai = 2:12:3, As = 2:12:22:3, Ae = x^x^ and Bi = -1102857 
/4003439a;ia;2 2:3 - 5464251/40034392:ia;2a;| + 2563669/40034392:^ 2:| + x^xl, B2 = -9223081 
/256164532:*a;^ - 18326919/25616453a:i2:^2:^ + 1933547/25616453a;*2^ + x^x^., 

B3 = -2617184886847/155355798195532:12:22:3 - 12918394932706/15535579819553xiX22:i^ + 2:^2:^, 

Bi = -26028972147097/1614908479876812:^2;^ - 135461875840584 

/16149084798768l2:j2:^2:^ + 2:^2:^, B5 = -2333331/36363072:^^2:22:^ - 1302976 

/36363072:i2:^2:3 + 2:12:22:3, Be = -11582471/3741935l2:i2:3 - 12629854 

/S7AW^blxlxlx3 - 4402342/124731172:52:^ + xixlxl, Bj = -xlx2X^ + xix^xf, 

Bs = -2:12:22:3 + xlx2xl, Bg = -2:12:3 + x\xlxl, Bio = -17362252580967/20864440421917x5*2:3 

- 3502187840950/208644404219172:!^ 2:22:3 + 2:'^2:|^, Bn = -x\x3 + xlxlx^. 

Figure 1: Motzkin's polynomial M (Eq. [T|) as Q2/Qi- 



This system was concocted by choosing Pi , P2 , ^3 somewhat haphazardly 
and then P4 = — (Pi+(3+(a;+52/)^)P2+P3+l+a;^), which guaranteed the system 
had no solution. The initial 130 constraints yield a search space of dimension 
145, and after four round of numeric solving one gets an unsatisfiability witness 
(sums of squares Qj such that X]j=i PjQj + *95 =0). Total computation time 
was 4.4 s. Even though there existed a simple solution (note the above formula 
for P4), our algorithm provided a lengthy one, with large coefficients (and thus 
unfit for inclusion here). 

Motzkin's polynomial M (Eq. [T|) cannot be expressed as a sum of squares, 
but it can be expressed as a quotient of two sums of squares. We solved 
M.Qi — Q2 = for sums of squares Qi and Q2 built from homogeneous mono- 
mials of respective total degrees 3 and 6 — lesser degrees yield no solutions 
(Fig. [1} . The equality relation over the polynomials yields 66 constraints over 
the matrix coefficients and a search space of dimension 186. Four cycles of SDP 
programming and LLL are then needed, total computation time was 4.1 s. 

We exhibited witnesses that each of the 8 semidefinite positive forms listed 
by [23| . which are not sums of squares of polynomials, are quotients of sums of 
squares (Motzkin's M, Robinson's R and /, Choi and Lam's F, Q, S, H and 
Schmiidgen's q). These examples include polynomials with up to 6 variables 
and search spaces up to dimension 1155. We did likewise with delzell, laxlax 
and leepstarr2 from 1^. The maximal computation time was 7'. 

We then converted these witnesses into Coq proofs of nonnegativity using a 
simple Sage script. These proofs use the Ring tactic, which checks for polyno- 
mial identity. Most proofs run within a few seconds, though laxlax takes 7'39" 
and Robinson's / 5'07"; the witness for leepstarr2 is too large for the parser. 
We also exhibited a witness that the Vorl polynomial cited by 24| is a sum of 
squares. 

John Harrison kindly provided us with a collection of 14 problems that his 
system [llj could not find witnesses for. These problems generally have the 
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form Pi>OA---AP„>0^i?>0. In order to prove such implication, we 
looked for witnesses consisting of sums of squares {Qi, . . . , Qm '5i?),such that 

QjPj + QrR — with Qji ^ 0, and thus R = — ^ . In some cases, it 
was necessary to use the products Jli -^T* ^o^' ^ {0,1}" instead of the Pi. 
We could find witnesses for all those problems[lll- though for some of them, the 
witnesses are very large, taking up megabytes. Since these searches were done 
without making use of symmetries in the problem, it is possible that more clever 
techniques could find smaller witnesses. 

5 Conclusion and further works 

We have described a method for solving SDP problems in rational arithmetic. 
This method can be used to solve sums-of-squares problems even in geometri- 
cally degenerate cases. We illustrated this method with applications to proving 
the nonnegativity of polynomials, or the unsatisfiability of systems of polyno- 
mial (in) equalities. The method then provides easily checkable proof witnesses, 
in the sense that checking the witness only entails performing polynomial arith- 
metic and applying a few simple mathematical lemmas. We have implemented 
the conversion of nonnegativeness witnesses to Coq proofs. A more ambitious 
implementation, mapping Coq real arithmetic proofs goals to Positivstellensatz 
problems through the Psatz tactic from the MicroMega package then map- 
ping Positivstellensatz witnesses back to proofs, is underway. 

One weakness of the method is that it tends to provide "unnatural" witnesses 
— they tend to have very large coefficients. These are machine-checkable but 
provide little insights to the reader. An alternative would be to provide the 
matrices and some additional data (such as their minimal polynomial) and have 
the checker verify that they are semidefinite positive; but this requires formally 
proving, once and for all, some non-trivial results on polynomials, symmetric 
matrices and eigenvalues (e.g. the Cayley-Hamilton theorem), as well as possibly 
performing costly computations, e.g. evaluating a matrix polynomial. 

A more serious limitation for proofs of unsatisfiability is the very high cost 
of application of the Positivstellensatz. There is the exponential number of 
polynomials to consider, and the unknown number of monomials. It would 
be very interesting if there could be some simple results, similar to the Newton 
polytope approach, for reducing the dimension of the search space or the number 
of polynomials to consider. Another question is whether it is possible to define 
SDP problems from Positivstellensatz equations for which the spectrahedron 
has rational points only at its relative boundary. 

While our method performed well on examples, and is guaranteed to provide 
a correct answer if it provides one, we have supplied no completeness proof — 
that is, we have not proved that it necessarily provides a solution if there is one. 
This is due to the use of fioating-point computations. One appreciable result 
would be that a solution should be found under the assumption that floating- 

^^A 7z archive is given at |http : //bit ■ ly/hM7HW3] 
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point computations are precise up to e, for a value of e and the various scaling 
factors in the algorithm depending on the values in the problem or the solution. 

It seems possible to combine our reduction method based on LLL with the 
Newton iterations suggested by 12, 1^, as an improvement over their strategy 
for detection of useless monomials and reduction of the search space. Again, 
further experiment is needed. 
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A Gaussian Reduction and Positive Semidefi- 
niteness 

An algorithm for transforming a semidefinite positive matrix into a "sum of 
squares" form, also known as Gaussian reduction: 

for i :=1 to n do 
begin 

i f m[ i , i ] < then 

throw non_positive_semidefinite 
i f m[ i , i ] =0 then 

i f m. row ( i ) <> 
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throw non_positive_semidefinite 

else 
begin 

V m. row ( i ) / m[ i , i ] 
output . append (ni[ i , i], v) 
m:=m — m[i, i] *v.transpose() *v 
end 

end 

Suppose that the entrance of iteration i, m[i . . .n,i . . .n] is positive semidef- 
inite. If mj,j = 0, then the ith base vector is in the isotropic cone of the matrix, 
thus of its kernel, and the row i must be zero. Otherwise, TOj ^ > 0. By adding e 
to the diagonal of the matrix, we would have a positive definite matrix and thus 
the output of the loop iteration would also be positive definite, as above. By 
e — > and the fact that the set of positive semidcfinitc matrices is topologically 
closed, then the output of the loop iteration is also positive semidefinite. 

The output variable is then a list of couples (cj, Vi) such that > and the 
original matrix m is equal to Civfvi (with Vi row vectors). Otherwise said, 
for any row vector u, umuF = Vi)"^ . 
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